********************************************************************************
* LPDENSITY STATA PACKAGE -- lpdensity mata functions
* Authors: Matias D. Cattaneo, Michael Jansson, Xinwei Ma
********************************************************************************
*!version 2.4 2023-01-21

** do lpdensity_fn.ado

capture mata mata drop lpdensity_fn()

mata
real matrix lpdensity_fn(real colvector data,
                         real colvector grid,
						 real colvector bw,
						 real scalar p, real scalar q, real scalar v,
						 string scalar kernel,
						 real colvector Cweights,
						 real colvector Pweights,
						 real scalar massPoints,
						 real scalar scale,
						 real scalar size,
						 real scalar CIuniform,
						 real scalar CIsimul){
					 
ii = order(data, 1)
data = data[ii, 1]
Cweights = Cweights[ii, 1]
Pweights = Pweights[ii, 1]
n  = length(data)
ng = length(grid)

dataUnique  = lpdensity_unique(data)
freqUnique  = dataUnique[., 2]
indexUnique = dataUnique[., 3]
dataUnique  = dataUnique[., 1]
nUnique     = length(dataUnique)

if (massPoints) {
	Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique)
} 
else{
	Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights)	
}

weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n
Cweights       = Cweights :/ sum(Cweights) :* n
Pweights       = Pweights :/ sum(Pweights) :* n

weights_normalUnique = runningsum(weights_normal)[indexUnique]
if (nUnique > 1) { 
	weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) 
}

CweightsUnique = runningsum(Cweights)[indexUnique]
if (nUnique > 1) { 
	CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) 
}

PweightsUnique = runningsum(Pweights)[indexUnique]
if (nUnique > 1) { 
	PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) 
}

Result = J(ng, 10, .)
iff_p = J(n, ng, .) // store influence functions
iff_q = iff_p

// loop over evaluation points
for (j=1; j<=ng; j++) {
    Result[j, 10] = j
	Result[j, 1] = grid[j]
	Result[j, 2] = bw[j]
	
	index_temp = (abs(data :- grid[j]) :<= bw[j])
	Result[j, 3] = sum(index_temp)
	
	// LHS and RHS variables
	if (massPoints) {
		Y_temp    = Fn[indexUnique]
		Xh_temp   = (data[indexUnique] :- grid[j]) :/ bw[j]

		// weights
		if (kernel == "triangular") {
			Kh_temp = ((1 :- abs(Xh_temp)) :/ bw[j]) :* index_temp[indexUnique]
		} 
		if (kernel == "uniform") {
			Kh_temp = (0.5 :/ bw[j]) :* index_temp[indexUnique]
		} 
		if (kernel == "epanechnikov") {
			Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ bw[j]) :* index_temp[indexUnique]
		}
	
		Xh_p_temp = J(nUnique, p+1, .)
		Xh_p_Kh_Pweights_temp = Xh_p_temp
	
		for (jj=1; jj<=p+1; jj++) {
			Xh_p_temp[., jj] = Xh_temp:^(jj-1)
			Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique
		}
	
	
		XhKhXh_inv = invsym(cross(Xh_p_temp, Xh_p_Kh_Pweights_temp) :/ n)
	
		// point estimate
		Result[j, 4] = factorial(v) :* (XhKhXh_inv * cross(Xh_p_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n
	
		// standard error estimate
		F_Xh_p_Kh_temp = cross(Y_temp, Xh_p_Kh_Pweights_temp) :/ n
	
		G = J(n, p - 0 + 1, .)
		for (jj=1; jj<=p-0+1; jj++) {
			G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal
		}
	
		iff_p[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))'
	}
	else {
		Y_temp    = Fn
		Xh_temp   = (data :- grid[j]) :/ bw[j]

		// weights
		if (kernel == "triangular") {
			Kh_temp = ((1 :- abs(Xh_temp)) :/ bw[j]) :* index_temp
		} 
		if (kernel == "uniform") {
			Kh_temp = (0.5 :/ bw[j]) :* index_temp
		} 
		if (kernel == "epanechnikov") {
			Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ bw[j]) :* index_temp
		}
	
		Xh_p_temp = J(n, p+1, .)
		Xh_p_Kh_Pweights_temp = Xh_p_temp
	
		for (jj=1; jj<=p+1; jj++) {
			Xh_p_temp[., jj] = Xh_temp:^(jj-1)
			Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights
		}
	
	
		XhKhXh_inv = invsym(cross(Xh_p_temp, Xh_p_Kh_Pweights_temp) :/ n)
	
		// point estimate
		Result[j, 4] = factorial(v) :* (XhKhXh_inv * cross(Xh_p_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n
	
		// standard error estimate
		F_Xh_p_Kh_temp = cross(Y_temp, Xh_p_Kh_Pweights_temp) :/ n
	
		G = J(n, p - 0 + 1, .)
		for (jj=1; jj<=p-0+1; jj++) {
			G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal
		}
	
		iff_p[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))'
	}
	
		
	// robust CI
	if (q > p) {
		if (massPoints) {
			Xh_q_temp = J(nUnique, q+1, .)
			Xh_q_Kh_Pweights_temp = Xh_q_temp
	
			for (jj=1; jj<=q+1; jj++) {
				Xh_q_temp[., jj] = Xh_temp:^(jj-1)
				Xh_q_Kh_Pweights_temp[., jj] = Xh_q_temp[., jj] :* Kh_temp :* PweightsUnique
			}
			XhKhXh_inv = invsym(cross(Xh_q_temp, Xh_q_Kh_Pweights_temp) :/ n)
		
			// point estimate
			Result[j, 5] = factorial(v) :* (XhKhXh_inv * cross(Xh_q_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n
	
			// standard error estimate
			F_Xh_q_Kh_temp = cross(Y_temp, Xh_q_Kh_Pweights_temp) :/ n
		
			G = J(n, q - 0 + 1, .)
			for (jj=1; jj<=q-0+1; jj++) {
				G[., jj] = (lpdensity_rep((runningsum(Xh_q_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_q_Kh_temp[jj]) :* weights_normal
			}
			iff_q[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))'
		}
		else {
			Xh_q_temp = J(n, q+1, .)
			Xh_q_Kh_Pweights_temp = Xh_q_temp
	
			for (jj=1; jj<=q+1; jj++) {
				Xh_q_temp[., jj] = Xh_temp:^(jj-1)
				Xh_q_Kh_Pweights_temp[., jj] = Xh_q_temp[., jj] :* Kh_temp :* Pweights
			}
			XhKhXh_inv = invsym(cross(Xh_q_temp, Xh_q_Kh_Pweights_temp) :/ n)
		
			// point estimate
			Result[j, 5] = factorial(v) :* (XhKhXh_inv * cross(Xh_q_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n
	
			// standard error estimate
			F_Xh_q_Kh_temp = cross(Y_temp, Xh_q_Kh_Pweights_temp) :/ n
		
			G = J(n, q - 0 + 1, .)
			for (jj=1; jj<=q-0+1; jj++) {
				G[., jj] = ((runningsum(Xh_q_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_q_Kh_temp[jj]) :* weights_normal
			}
			iff_q[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))'
		}
	}
}

CovMat_p = cross(iff_p, iff_p) :/ n
CovMat_q = cross(iff_q, iff_q) :/ n

Result[., 6] = sqrt(abs(diagonal(CovMat_p)))
Result[., 7] = sqrt(abs(diagonal(CovMat_q)))
	
Result[., 4..7] = Result[., 4..7] :* scale

if (CIuniform == 0) {
	critVal = invnormal(1-size/2)
}
else {
	if (q > p) {
		CovMat_q = CovMat_q :* scale^2
		for (j1=1; j1 <= ng; j1++) {
			for (j2=1; j2 <= ng; j2++) {
				CovMat_q[j1, j2] = CovMat_q[j1, j2] / Result[j1, 7] / Result[j2, 7]
			}
		}
		normRVs = cholesky(CovMat_q) * rnormal(ng, CIsimul, 0, 1)
		critVal = lpdensity_quantile(colmax(abs(normRVs))', 1-size)
	}
	else {
		CovMat_p = CovMat_p :* scale^2
		for (j1=1; j1 <= ng; j1++) {
			for (j2=1; j2 <= ng; j2++) {
				CovMat_p[j1, j2] = CovMat_p[j1, j2] / Result[j1, 6] / Result[j2, 6]
			}
		}
		normRVs = cholesky(CovMat_p) * rnormal(ng, CIsimul, 0, 1)
		critVal = lpdensity_quantile(colmax(abs(normRVs))', 1-size)
	}
}
	
if (q > p) {
	Result[., 8] = Result[., 5] - critVal * Result[., 7]
	Result[., 9] = Result[., 5] + critVal * Result[., 7]
}
if (q == p) {
	Result[., 8] = Result[., 4] - critVal * Result[., 6]
	Result[., 9] = Result[., 4] + critVal * Result[., 6]
}

return(Result)
	
}
mata mosave lpdensity_fn(), replace
end

********************************************************************************
* Extracting unique elements
********************************************************************************
capture mata mata drop lpdensity_unique()

mata
real matrix lpdensity_unique(real colvector x){
// Note: x should be in ascending order

n = length(x)

// x has one or no element
if (n == 0) {
	out = J(0, 3, .)
	return(out)
}
if (n == 1) {
	out = (x, 1, 1)
	return(out)
}

// else
numIndexTemp = selectindex(x[2..n] :!= x[1..(n-1)])
if (length(numIndexTemp) == 0) {
	numIndexLast  = J(1, 1, n)
	numIndexFirst = J(1, 1, 1)
} else {
	numIndexLast  = (numIndexTemp \ n)
	numIndexFirst = (1 \ numIndexTemp:+1)
}

freq = numIndexLast - numIndexFirst :+ 1
unique = x[numIndexLast]

out = (unique, freq, numIndexLast)
return(out)
}
/* // testing
mata: x = (1\ 1)
mata: lpdensity_unique(x)

mata: x = (1\ 2\ 2)
mata: lpdensity_unique(x)

mata: x = (1..100)'
mata: lpdensity_unique(x)

mata: x = (1\ 1\ 2\ 3\ 4\ 4\ 5\ 6\ 6\ 6\ 7)
mata: lpdensity_unique(x)

mata: x = (1\ 1\ 2\ 3\ 4\ 4\ 5\ 6\ 6\ 6\ 7\ 7)
mata: lpdensity_unique(x)

mata: x = (5)
mata: lpdensity_unique(x)
*/ 
mata mosave lpdensity_unique(), replace
end

********************************************************************************
* Select the minimum element
********************************************************************************
capture mata mata drop lpdensity_whichmin()

mata
real scalar lpdensity_whichmin(real colvector x){
	
if (length(x) == 0) {
	return(0)
}
else {
	return (order(x, 1)[1])	
}


}
/* // testing
mata: x = (1\ 1)
mata: lpdensity_whichmin(x)

mata: x = (1\ 2\ 2)
mata: lpdensity_whichmin(x)

mata: x = (3\ 2\ 2)
mata: lpdensity_whichmin(x)

mata: x = (3\ 2\ 3)
mata: lpdensity_whichmin(x)

mata: x = (3\ 2\ 1)
mata: lpdensity_whichmin(x)
*/ 
mata mosave lpdensity_whichmin(), replace
end

********************************************************************************
* Replicating each element in x y times (elementwise)
********************************************************************************
capture mata mata drop lpdensity_rep()

mata
real matrix lpdensity_rep(real colvector x, real colvector y){
// Warning: x and y should have the same length
//          y should contain strictly positive integers
//          y cannot be empty

nout = sum(y)
if (length(y) == 1) {
	out = J(nout, 1, x)
}
else {
	if (all(y :== 1)) {
		out = x
	}
	else {
		out = J(nout, 1, .)
	    out[1..y[1], 1] = J(y[1], 1, x[1])
	    indextemp = y[1]
	    for (i=2; i<=length(y); i++) {
			out[(indextemp+1)..(indextemp+y[i])] = J(y[i], 1, x[i])
			indextemp = indextemp + y[i]

		}
	}
}

return(out)
}
/* // testing
mata: lpdensity_rep(5, 1)
mata: lpdensity_rep(5, 3)
mata: lpdensity_rep((5\ 6), (1\ 1))
mata: lpdensity_rep((5\ 6), (1\ 2))
mata: lpdensity_rep((5\ 6), (2\ 1))
mata: lpdensity_rep((5\ 6), (2\ 3))
mata: lpdensity_rep((5\ 6\ 7\ 8), (1\ 1\ 1\ 1))
mata: lpdensity_rep((5\ 6\ 7\ 8), (1\ 2\ 2\ 3))
mata: lpdensity_rep((5\ 6\ 7\ 8), (3\ 1\ 2\ 1))
mata: lpdensity_rep((5\ 6\ 7\ 8), (3\ 2\ 2\ 3))
*/ 
mata mosave lpdensity_rep(), replace
end

********************************************************************************
* Normal DGPs
********************************************************************************
capture mata mata drop lpdensity_norm_dgps()

mata
real scalar lpdensity_normDgp(real scalar y, real scalar v, real scalar mu, real scalar sigma){

x = (y - mu) / sigma // normalized evaluation point

if (v == 0 ) out = normal(x) // the cdf

if (v > 0) {
	// compute hermite polynomials
	if (v == 1) H = 1
	if (v == 2) H = x
	if (v > 2) {
		temp = J(v, 1, .)
		temp[1] = 1; temp[2] = x
		for (j=3; j<=v; j++) temp[j] = x * temp[j-1] - (j-2) * temp[j-2]
		H = temp[v]
	}
	// get the sign correct
	H = H * (-1)^(v + 1)	
	// get the scaling correct
	out = H * normalden(x) / sigma^v
}

return(out)

}
mata mosave lpdensity_normDgp(), replace
end

********************************************************************************
* empirical quantile
********************************************************************************
capture mata mata drop lpdensity_quantile()

mata
real scalar lpdensity_quantile(real colvector x, real scalar p){

x = sort(x, 1)
n = length(x)
q = ceil(p * n)
if (q < 1) q = 1
if (q > n) q = n
out = x[q]

return(out)

}
mata mosave lpdensity_quantile(), replace
end

********************************************************************************
* S matrix generation
********************************************************************************
capture mata mata drop lpdensity_Sgen()

mata
real matrix lpdensity_Sgen(real scalar p, string scalar kernel){

if (kernel == "uniform") {
	out =  (1,0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909\ 
			0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0\ 
			0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769\ 
			0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0\ 
			0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667\ 
			0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0\ 
			0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647\ 
			0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0\ 
			0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684\ 
			0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0\ 
			0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476) 
}
if (kernel == "triangular") {
	out =  (1,0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831\ 
			0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0\ 
			0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887\ 
			0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0\ 
			0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318\ 
			0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0\ 
			0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977\ 
			0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0\ 
			0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991\ 
			0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0\ 
			0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433)
} 
if (kernel == "epanechnikov") {
	out =  (1,0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021\ 
			0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0\ 
			0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154\ 
			0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0\ 
			0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529\ 
			0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0\ 
			0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443\ 
			0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0\ 
			0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812\ 
			0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0\ 
			0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236)
}
return(out[1..(p+1), 1..(p+1)])

}
mata mosave lpdensity_Sgen(), replace
end

********************************************************************************
* T matrix generation
********************************************************************************
capture mata mata drop lpdensity_Tgen()

mata
real matrix lpdensity_Tgen(real scalar p, string scalar kernel){

if (kernel == "uniform") {
	out =  (0.5,0,0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455\ 
			0,0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0\ 
			0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385\ 
			0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0\ 
			0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333\ 
			0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0\ 
			0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824\ 
			0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0\ 
			0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842\ 
			0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842,0\ 
			0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842,0,0.0238095238095238) 
}
if (kernel == "triangular") {
	out =  (0.666666666666667,0,0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814\ 
			0,0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0\ 
			0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681\ 
			0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0\ 
			0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244\ 
			0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0\ 
			0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965\ 
			0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0\ 
			0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208\ 
			0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208,0\ 
			0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208,0,0.000376435159043855)
} 
if (kernel == "epanechnikov") {
	out =  (0.6,0,0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042\ 
			0,0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0\ 
			0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683\ 
			0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0\ 
			0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889\ 
			0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0\ 
			0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492\ 
			0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0\ 
			0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932\ 
			0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932,0\ 
			0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932,0,0.000745341614906832)
}
return(out[1..(p+1), 1..(p+1)])

}
mata mosave lpdensity_Tgen(), replace
end

********************************************************************************
* C matrix generation
********************************************************************************
capture mata mata drop lpdensity_Cgen()

mata
real matrix lpdensity_Cgen(real scalar k, real scalar p, string scalar kernel){

if (kernel == "uniform") {
	out =  (1,0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0\ 
			0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647\ 
			0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0\ 
			0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684\ 
			0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0\ 
			0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476\ 
			0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0\ 
			0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652\ 
			0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0\ 
			0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0,0.04\ 
			0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0,0.04,0) 
}
if (kernel == "triangular") {
	out =  (1,0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0\ 
			0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977\ 
			0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0\ 
			0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991\ 
			0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0\ 
			0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433\ 
			0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0\ 
			0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971\ 
			0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0\ 
			0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0,0.00307692307692308\ 
			0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0,0.00307692307692308,0)
} 
if (kernel == "epanechnikov") {
	out =  (1,0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0\ 
			0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443\ 
			0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0\ 
			0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812\ 
			0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0\ 
			0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236\ 
			0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0\ 
			0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783\ 
			0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0\ 
			0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0,0.00444444444444445\ 
			0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0,0.00444444444444445,0)
}
return(out[1..(p+1), k+1])

}
mata mosave lpdensity_Cgen(), replace
end

********************************************************************************
* G matrix generation
********************************************************************************
capture mata mata drop lpdensity_Ggen()

mata
real matrix lpdensity_Ggen(real scalar p, string scalar kernel){

if (kernel == "uniform") {
	out =  (-0.333333333333333,0.166666666666667,-0.133333333333333,0.1,-0.0857142857142857,0.0714285714285714,-0.0634920634920635,0.0555555555555556,-0.0505050505050505,0.0454545454545455,-0.041958041958042\ 
			0.166666666666667,0.0666666666666667,0.0555555555555556,0.0380952380952381,0.0333333333333333,0.0264550264550265,0.0238095238095238,0.0202020202020202,0.0185185185185185,0.0163170163170163,0.0151515151515152\ 
			-0.133333333333333,0.0555555555555556,-0.0476190476190476,0.0333333333333333,-0.0296296296296296,0.0238095238095238,-0.0216450216450216,0.0185185185185185,-0.0170940170940171,0.0151515151515152,-0.0141414141414141\ 
			0.1,0.0380952380952381,0.0333333333333333,0.0222222222222222,0.02,0.0155844155844156,0.0142857142857143,0.011965811965812,0.0111111111111111,0.0096969696969697,0.00909090909090909\ 
			-0.0857142857142857,0.0333333333333333,-0.0296296296296296,0.02,-0.0181818181818182,0.0142857142857143,-0.0131868131868132,0.0111111111111111,-0.0103703703703704,0.00909090909090909,-0.00855614973262032\ 
			0.0714285714285714,0.0264550264550265,0.0238095238095238,0.0155844155844156,0.0142857142857143,0.010989010989011,0.0102040816326531,0.00846560846560847,0.00793650793650794,0.00687547746371276,0.0064935064935065\ 
			-0.0634920634920635,0.0238095238095238,-0.0216450216450217,0.0142857142857143,-0.0131868131868132,0.0102040816326531,-0.00952380952380953,0.00793650793650794,-0.00746965452847806,0.0064935064935065,-0.00615174299384826\ 
			0.0555555555555556,0.0202020202020202,0.0185185185185185,0.011965811965812,0.0111111111111111,0.00846560846560847,0.00793650793650794,0.0065359477124183,0.00617283950617284,0.00531632110579479,0.00505050505050505\ 
			-0.0505050505050505,0.0185185185185185,-0.0170940170940171,0.0111111111111111,-0.0103703703703704,0.00793650793650794,-0.00746965452847806,0.00617283950617284,-0.00584795321637427,0.00505050505050505,-0.00481000481000481\ 
			0.0454545454545455,0.0163170163170163,0.0151515151515152,0.0096969696969697,0.00909090909090909,0.00687547746371276,0.00649350649350649,0.00531632110579479,0.00505050505050505,0.00432900432900433,0.00413223140495868\ 
			-0.041958041958042,0.0151515151515152,-0.0141414141414141,0.00909090909090909,-0.00855614973262032,0.00649350649350649,-0.00615174299384826,0.00505050505050505,-0.00481000481000481,0.00413223140495868,-0.00395256916996048) 
}
if (kernel == "triangular") {
	out =  (-0.233332866779462,0.0833330508905063,-0.0531743640785073,0.0333333960560459,-0.0243382259715177,0.0178568590991684,-0.0140329271889015,0.0111109568092491,-0.00914256712672596,0.00757561833195314,-0.00643232469904509\ 
			0.0833333430031619,0.0206349536543496,0.0138889822607069,0.00747358076709329,0.00555558002800799,0.00376385059786869,0.00297624095083857,0.00224683629513331,0.00185187157177693,0.00148741419069681,0.00126263242255087\ 
			-0.0531747120633688,0.0138889436189969,-0.00992067054772658,0.00555558570938115,-0.00427612357529663,0.00297621392905644,-0.00240733368919841,0.00185187029011056,-0.00155068139981584,0.00126263712216438,-0.00108438068750666\ 
			0.0333333291761154,0.00747354338303244,0.0055555555888295,0.00282828659007681,0.00222222610846186,0.00145548770075321,0.00119047711034882,0.000879119696771778,0.000740741774432789,0.000586005979164481,0.00050505152540861\ 
			-0.0243386235494463,0.00555555541211971,-0.00427609423962765,0.0022222251541454,-0.00178710387992154,0.00119047600313231,-0.000989883016025128,0.000740740651962696,-0.00063180824060111,0.000505050436046404,-0.000439378283083419\ 
			0.0178571428772336,0.0037638288080394,0.00297619046584015,0.00145549115212934,0.00119047780581393,0.000758765013497074,0.000637755044256656,0.000461990375582904,0.000396825380503709,0.000309560314164898,0.000270562752861001\ 
			-0.014033189060362,0.00297619047992443,-0.0024073149074495,0.00119047782568061,-0.000989884403822643,0.000637755064921552,-0.00054271705691498,0.000396825400688911,-0.000344133756729751,0.000270562772022687,-0.000238285106683998\ 
			0.01111111111111,0.00224682724469665,0.00185185185231585,0.000879122147763556,0.000740741756863993,0.000461990319945608,0.000396825372642262,0.000282842182321384,0.000246913581918653,0.000190248347383829,0.000168350168854074\ 
			-0.00914270914122569,0.00185185185168499,-0.00155067155074166,0.000740741756146379,-0.000631809109714714,0.000396825371926143,-0.000344133734564222,0.000246913581232753,-0.00021720969177002,0.000168350168212254,-0.000149908647886651\ 
			0.00757575757574178,0.00148740148756518,0.00126256257569384,0.000586007961071109,0.000505051197433514,0.000309560273844802,0.000270562753630222,0.000190248349061672,0.00016835016904926,0.000128330167831491,0.000114784205611466\ 
			-0.00643245643257991,0.00126256257573434,-0.00108431925428119,0.000505051197475794,-0.000439378889047697,0.000270562753671301,-0.000238285092406906,0.000168350169087952,-0.000149908648550576,0.000114784205647198,-0.000103205972768449)
} 
if (kernel == "epanechnikov") {
	out =  (-0.257142857142857,0.1,-0.0666666666666667,0.0428571428571429,-0.032034632034632,0.0238095238095238,-0.018981018981019,0.0151515151515152,-0.0125874125874126,0.0104895104895105,-0.00896750308515014\ 
			0.1,0.0285714285714286,0.02,0.0112554112554113,0.00857142857142857,0.00592740592740593,0.00476190476190476,0.00363636363636364,0.00303030303030303,0.00245166598107775,0.0020979020979021\ 
			-0.0666666666666667,0.02,-0.0147186147186147,0.00857142857142857,-0.00672660672660673,0.00476190476190476,-0.0039027639027639,0.00303030303030303,-0.00256136020841903,0.0020979020979021,-0.00181428478642101\ 
			0.0428571428571429,0.0112554112554113,0.00857142857142857,0.00459540459540459,0.0036734693877551,0.00246420246420246,0.00204081632653061,0.00152710035062976,0.0012987012987013,0.00103586815661119,0.000899100899100899\ 
			-0.032034632034632,0.00857142857142857,-0.00672660672660672,0.0036734693877551,-0.002997002997003,0.00204081632653061,-0.00171514759750054,0.0012987012987013,-0.00111669548202056,0.000899100899100899,-0.000787094374002821\ 
			0.0238095238095238,0.00592740592740593,0.00476190476190476,0.00246420246420246,0.00204081632653061,0.00133591898297781,0.00113378684807256,0.000833634889362443,0.000721500721500722,0.000568059391588803,0.0004995004995005\ 
			-0.018981018981019,0.00476190476190476,-0.0039027639027639,0.00204081632653061,-0.00171514759750054,0.00113378684807256,-0.000973020787262273,0.000721500721500722,-0.000629917038585769,0.0004995004995005,-0.000442294982994558\ 
			0.0151515151515152,0.00363636363636364,0.00303030303030303,0.00152710035062976,0.0012987012987013,0.000833634889362443,0.000721500721500721,0.000522697117124362,0.000459136822773186,0.000357384796744064,0.000317863954227591\ 
			-0.0125874125874126,0.00303030303030303,-0.00256136020841903,0.0012987012987013,-0.00111669548202056,0.000721500721500721,-0.000629917038585769,0.000459136822773186,-0.000406153724231527,0.000317863954227591,-0.000284353327831589\ 
			0.0104895104895105,0.00245166598107775,0.0020979020979021,0.00103586815661119,0.000899100899100899,0.000568059391588803,0.000499500499500499,0.000357384796744064,0.000317863954227591,0.000244972418885462,0.000220059660619101\ 
			-0.00896750308515014,0.0020979020979021,-0.00181428478642101,0.000899100899100899,-0.000787094374002822,0.000499500499500499,-0.000442294982994558,0.000317863954227591,-0.000284353327831589,0.000220059660619101,-0.000198641937772373)
}
return(out[1..(p+1), 1..(p+1)])

}
mata mosave lpdensity_Ggen(), replace
end

********************************************************************************
* mse-rot bandwidth
********************************************************************************
capture mata mata drop lpdensity_bwROT()

mata
real matrix lpdensity_bwROT(real colvector data,
							real colvector grid,
							real scalar p, real scalar v,
							string scalar kernel,
							real colvector Cweights,
							real colvector Pweights,
							real scalar massPoints,
							real scalar stdVar,
							real scalar regularize,
							real scalar nLocalMin,
							real scalar nUniqueMin){
							
data = sort(data, 1)
					
n 	= length(data)
ng 	= length(grid)

dataUnique = lpdensity_unique(data)[., 1]
nUnique = length(dataUnique)

if (stdVar) {
	center_temp = mean(data)
	scale_temp  = sqrt(variance(data))
	
	data = (data :- center_temp) :/ scale_temp
	dataUnique = (dataUnique :- center_temp) :/ scale_temp
	grid = (grid :- center_temp) :/ scale_temp
}

// estimate a normal reference model
mean_hat = sum(data :* Cweights :* Pweights) / sum(Cweights :* Pweights)
sd_hat   = sqrt(sum(Cweights :* Pweights :* (data :- mean_hat):^2) / sum(Cweights :* Pweights))

// bias estimate, no rate added
bias_dgp = J(ng, 2, .)
for (j=1; j<=ng; j++) {
	bias_dgp[j, 1] = lpdensity_normDgp(grid[j], p+1, mean_hat, sd_hat) / factorial(p+1) * factorial(v)
	bias_dgp[j, 2] = lpdensity_normDgp(grid[j], p+2, mean_hat, sd_hat) / factorial(p+2) * factorial(v) +
		bias_dgp[j, 1] * lpdensity_normDgp(grid[j], 2, mean_hat, sd_hat) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat)
 }
 
S  = lpdensity_Sgen(     p, kernel=kernel)
C1 = lpdensity_Cgen(p+1, p, kernel=kernel)
C2 = lpdensity_Cgen(p+2, p, kernel=kernel)
G  = lpdensity_Ggen(     p, kernel=kernel)
S2 = lpdensity_Tgen(     p, kernel=kernel)

bias_dgp[., 1] = bias_dgp[., 1] :* (invsym(S) * C1)[v+1, .]
bias_dgp[., 2] = bias_dgp[., 2] :* (invsym(S) * C2)[v+1, .]

// variance estimate, sample size added
sd_dgp = J(ng, 1, .)
if (v > 0) {
	for (j=1; j<=ng; j++) {
		sd_dgp[j, 1] = factorial(v) * sqrt(lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / n)
	}
	sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * G * invsym(S))[v+1, v+1]))
} 
else {
	for (j=1; j<=ng; j++) {
		sd_dgp[j, 1] = factorial(v) * sqrt(
		lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat) * (1 - lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat)) /
		lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / (0.5 * n^2))
	}
	sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * S2 * invsym(S))[v+1, v+1]))
}

// bandwidth
h = J(ng, 1, .)

for (j=1; j<=ng; j++) {

	S = optimize_init()
	optimize_init_evaluator(S, &lpdensity_optFunc())
	optimize_init_params(S, 1e-8)
	optimize_init_argument(S, 1, bias_dgp[j, .])
	optimize_init_argument(S, 2, sd_dgp[j, .])
	optimize_init_argument(S, 3, p)
	optimize_init_argument(S, 4, v)
	optimize_init_argument(S, 5, "mse-rot")
	optimize_init_which(S, "min")
	optimize_init_evaluatortype(S, "gf0")
	optimize_init_tracelevel(S, "none")
	optimize_init_conv_maxiter(S, 20)
	(void) optimize(S)
	h[j] = optimize_result_params(S)
	
	if (regularize == 1) {
		if (nLocalMin > 0) {
			h[j] = max((h[j], sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))]))
		}
		if (nUniqueMin > 0) {
			h[j] = max((h[j], sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))]))
		}
	}
	h[j] = min((  h[j], max(  abs(dataUnique :- grid[j])  )  ))
}

if (stdVar) {
	h = h :* scale_temp
}

return(h) 

}
mata mosave lpdensity_bwROT(), replace
end

********************************************************************************
* imse-rot bandwidth
********************************************************************************
capture mata mata drop lpdensity_bwIROT()

mata
real scalar lpdensity_bwIROT(real colvector data,
							 real colvector grid,
							 real scalar p, real scalar v,
							 string scalar kernel,
							 real colvector Cweights,
							 real colvector Pweights,
							 real scalar massPoints,
							 real scalar stdVar,
							 real scalar regularize,
							 real scalar nLocalMin,
							 real scalar nUniqueMin){
						
data = sort(data, 1)
					
n 	= length(data)
ng 	= length(grid)

dataUnique = lpdensity_unique(data)[., 1]
nUnique = length(dataUnique)

if (stdVar) {
	center_temp = mean(data)
	scale_temp  = sqrt(variance(data))

	data = (data :- center_temp) :/ scale_temp
	dataUnique = (dataUnique :- center_temp) :/ scale_temp
	grid = (grid :- center_temp) :/ scale_temp
}

// estimate a normal reference model
mean_hat = sum(data :* Cweights :* Pweights) / sum(Cweights :* Pweights)
sd_hat   = sqrt(sum(Cweights :* Pweights :* (data :- mean_hat):^2) / sum(Cweights :* Pweights))

// bias estimate, no rate added
bias_dgp = J(ng, 2, .)
for (j=1; j<=ng; j++) {
	bias_dgp[j, 1] = lpdensity_normDgp(grid[j], p+1, mean_hat, sd_hat) / factorial(p+1) * factorial(v)
	bias_dgp[j, 2] = lpdensity_normDgp(grid[j], p+2, mean_hat, sd_hat) / factorial(p+2) * factorial(v) +
		bias_dgp[j, 1] * lpdensity_normDgp(grid[j], 2, mean_hat, sd_hat) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat)
 }
 
S  = lpdensity_Sgen(     p, kernel=kernel)
C1 = lpdensity_Cgen(p+1, p, kernel=kernel)
C2 = lpdensity_Cgen(p+2, p, kernel=kernel)
G  = lpdensity_Ggen(     p, kernel=kernel)
S2 = lpdensity_Tgen(     p, kernel=kernel)

bias_dgp[., 1] = bias_dgp[., 1] :* (invsym(S) * C1)[v+1, .]
bias_dgp[., 2] = bias_dgp[., 2] :* (invsym(S) * C2)[v+1, .]

// variance estimate, sample size added
sd_dgp = J(ng, 1, .)
if (v > 0) {
	for (j=1; j<=ng; j++) {
		sd_dgp[j, 1] = factorial(v) * sqrt(lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / n)
	}
	sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * G * invsym(S))[v+1, v+1]))
} 
else {
	for (j=1; j<=ng; j++) {
		sd_dgp[j, 1] = factorial(v) * sqrt(
		lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat) * (1 - lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat)) /
		lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / (0.5 * n^2))
	}
	sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * S2 * invsym(S))[v+1, v+1]))
}

// bandwidth
S = optimize_init()
optimize_init_evaluator(S, &lpdensity_optFunc())
optimize_init_params(S, 1e-8)
optimize_init_argument(S, 1, bias_dgp)
optimize_init_argument(S, 2, sd_dgp)
optimize_init_argument(S, 3, p)
optimize_init_argument(S, 4, v)
optimize_init_argument(S, 5, "imse-rot")
optimize_init_which(S, "min")
optimize_init_evaluatortype(S, "gf0")
optimize_init_tracelevel(S, "none")
optimize_init_conv_maxiter(S, 20)
(void) optimize(S)
h = optimize_result_params(S)

if (regularize == 1) {
	for (j=1; j<=ng; j++) {
		if (nLocalMin > 0) {
			h = max((h, sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))]))
		}
		if (nUniqueMin > 0) {
			h = max((h, sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))]))
		}
	}
}

h = min(( h, max((  abs(max(dataUnique)-min(grid)), abs(min(dataUnique)-max(grid))  ))  ))

if (stdVar) {
	h = h:* scale_temp
}

return(h) 

}
mata mosave lpdensity_bwIROT(), replace
end

********************************************************************************
* mse-dpi bandwidth
********************************************************************************
capture mata mata drop lpdensity_bwMSE()

mata
real matrix lpdensity_bwMSE(real colvector data,
							real colvector grid,
							real scalar p, real scalar v,
							string scalar kernel,
							real colvector Cweights,
							real colvector Pweights,
							real scalar massPoints,
							real scalar stdVar,
							real scalar regularize,
							real scalar nLocalMin,
							real scalar nUniqueMin){

					
ii = order(data, 1)
data = data[ii, 1]
Cweights = Cweights[ii, 1]
Pweights = Pweights[ii, 1]					
n  = length(data)
ng = length(grid)

dataUnique  = lpdensity_unique(data)
freqUnique  = dataUnique[., 2]
indexUnique = dataUnique[., 3]
dataUnique  = dataUnique[., 1]
nUnique     = length(dataUnique)

if (stdVar) {
	center_temp = mean(data)
	scale_temp  = sqrt(variance(data))

	data = (data :- center_temp) :/ scale_temp
	dataUnique = (dataUnique :- center_temp) :/ scale_temp
	grid = (grid :- center_temp) :/ scale_temp
}

if (massPoints) {
	Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique)
} 
else{
	Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights)	
}

weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n
Cweights       = Cweights :/ sum(Cweights) :* n
Pweights       = Pweights :/ sum(Pweights) :* n

weights_normalUnique = runningsum(weights_normal)[indexUnique]
if (nUnique > 1) { 
	weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) 
}

CweightsUnique = runningsum(Cweights)[indexUnique]
if (nUnique > 1) { 
	CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) 
}

PweightsUnique = runningsum(Pweights)[indexUnique]
if (nUnique > 1) { 
	PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) 
}

// obtain preliminary bandwidth for estimating densities
// this is used for constructing preasymptotic matrices
h1  = lpdensity_bwIROT(data, grid, 2,   1,   kernel, Cweights, Pweights, 1, 1, 1, 20+2+1,   20+2+1)

// obtain preliminary bandwidth for estimating F_p+1
// this is used for constructing F_p+1
hp1 = lpdensity_bwIROT(data, grid, p+2, p+1, kernel, Cweights, Pweights, 1, 1, 1, 20+p+2+1, 20+p+2+1)

// obtain preliminary bandwidth for estimating F_p+2
// this is used for constructing F_p+2
hp2 = lpdensity_bwIROT(data, grid, p+3, p+2, kernel, Cweights, Pweights, 1, 1, 1, 20+p+3+1, 20+p+3+1)

dgp_hat = J(ng, 2, .) // Fp+1 and Fp+2 with normalization constants
const_hat = J(ng, 3, .)
h = J(ng, 1, .)

for (j=1; j<=ng; j++) {

	// estimate F_p+2
	index_temp = (abs(data :- grid[j]) :<= hp2)
	Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp2 
	Xh_p_temp = J(sum(index_temp), p+4, 1)	
	Xh_p_temp[., 2] = Xh_temp

	for (jj=3; jj<=p+4; jj++) { 
		Xh_p_temp[., jj] = Xh_temp :^ (jj-1)
	}
		
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp)) :/ hp2
	if (kernel == "uniform") 		Kh_temp = 0.5 / hp2
	if (kernel == "epanechnikov") 	Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp2
	
	Kh_temp = select(Pweights, index_temp) :* Kh_temp
	Y_temp = select(Fn, index_temp)
	dgp_hat[j, 2] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+3, 1]  / (hp2^(p+2))

	// estimate F_p+1
    index_temp = (abs(data :- grid[j]) :<= hp1)
	Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp1 
	Xh_p_temp = J(sum(index_temp), p+3, 1)
	Xh_p_temp[., 2] = Xh_temp
	
	for (jj=3; jj<=p+3; jj++) { 
		Xh_p_temp[., jj] = Xh_temp:^(jj-1)
	}
	
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp))  :/ hp1
	if (kernel == "uniform") 		Kh_temp = 0.5 / hp1
	if (kernel == "epanechnikov") 	Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp1
	
	Kh_temp = select(Pweights, index_temp) :* Kh_temp
	Y_temp = select(Fn, index_temp)
	dgp_hat[j, 1] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+2, 1] / (hp1^(p+1))

	// prepare for estimating matrices
    index_temp = (abs(data :- grid[j]) :<= h1)
    Xh_temp = (select(data, index_temp) :- grid[j]) :/ h1 
	
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp)) :/ h1
	if (kernel == "uniform") 		Kh_temp = 0.5 / h1
	if (kernel == "epanechnikov") 	Kh_temp = 0.75 :* (1 :- Xh_temp:^2) :/ h1
	Pweights_temp = select(Pweights, index_temp)

    // estimate Cp matrix
	Xh_p_temp = J(sum(index_temp), (2*p+1) - (p+1) + 1, 1)
	for (jj=1; jj<=(2*p+1)-(p+1)+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(p+jj)
	}
	C_p_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n

	// estimate Cp+1 matrix
	Xh_p_temp = J(sum(index_temp), (2*p+2) - (p+2) + 1, 1)
	for (jj=1; jj<=(2*p+2)-(p+2)+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(p+jj+1)
	}
	C_p1_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n
	
	// estimate S matirx
 	Xh_p_temp = J(sum(index_temp), p - 0 + 1, 1)
	for (jj=1; jj<=p-0+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(jj-1)
	}

	S_hat = cross(Xh_p_temp, Kh_temp :* Pweights_temp, Xh_p_temp) :/ n
	S_hat_inv = invsym(S_hat)
	
	// estimate G matrix
	if (v == 0) {
		G_hat = cross(Xh_p_temp, Kh_temp:^2 :* Pweights_temp, Xh_p_temp) :/ n
	} 
	else {
		if (massPoints) {
			Y_temp = Fn[indexUnique]
			Xh_temp = (dataUnique :- grid[j]) :/ h1 
			Xh_p_temp = J(nUnique, p - 0 + 1, 1)
			Xh_p_Kh_Pweights_temp = Xh_p_temp
		
			if (kernel == "triangular")  	Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp[indexUnique]
			if (kernel == "uniform") 		Kh_temp =  (0.5 / h1) :* index_temp[indexUnique]
			if (kernel == "epanechnikov") 	Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp[indexUnique]
		
			for (jj=1; jj<=p-0+1; jj++) {
				Xh_p_temp[., jj] = Xh_temp:^(jj-1)
				Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique
			}
		
			F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n
		
			G = J(n, p - 0 + 1, .)
		
			for (jj=1; jj<=p-0+1;jj++) {
				G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal
			}

			G_hat = cross(G, G) :/ n	
		}
		else {
			Y_temp = Fn
			Xh_temp = (data :- grid[j]) :/ h1 
			Xh_p_temp = J(n, p - 0 + 1, 1)
			Xh_p_Kh_Pweights_temp = Xh_p_temp
		
			if (kernel == "triangular")  	Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp
			if (kernel == "uniform") 		Kh_temp =  (0.5 / h1) :* index_temp
			if (kernel == "epanechnikov") 	Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp
		
			for (jj=1; jj<=p-0+1; jj++) {
				Xh_p_temp[., jj] = Xh_temp:^(jj-1)
				Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights
			}
		
			F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n
		
			G = J(n, p - 0 + 1, .)
		
			for (jj=1; jj<=p-0+1;jj++) {
				G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal
			}

			G_hat = cross(G, G) :/ n
		}
	}


	// now get all constants
    const_hat[j, 1] = factorial(v) * (S_hat_inv * C_p_hat )[v+1]
    const_hat[j, 2] = factorial(v) * (S_hat_inv * C_p1_hat)[v+1]

    if (v > 0) {
		const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (n * h1))
    } 
	else {
		temp_ii = min((  max((  mean(data :<= grid[j]), 1/n  )), 1 - 1/n ))
		const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (0.5 * n^2) * 
			h1 * temp_ii * (1 - temp_ii))
    }
	
	// now optimal bandwidth
	S = optimize_init()
	optimize_init_evaluator(S, &lpdensity_optFunc())
	optimize_init_params(S, 1e-8)
	optimize_init_argument(S, 1, dgp_hat[j, .])
	optimize_init_argument(S, 2, const_hat[j, .])
	optimize_init_argument(S, 3, p)
	optimize_init_argument(S, 4, v)
	optimize_init_argument(S, 5, "mse-dpi")
	optimize_init_which(S, "min")
	optimize_init_evaluatortype(S, "gf0")
	optimize_init_tracelevel(S, "none")
	optimize_init_conv_maxiter(S, 20)
	(void) optimize(S)
	h[j] = optimize_result_params(S)
	

	if (regularize == 1) {
		if (nLocalMin > 0) {
			h[j] = max((h[j], sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))]))
		}
		if (nUniqueMin > 0) {
			h[j] = max((h[j], sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))]))
		}
	}
	h[j] = min((  h[j], max(  abs(dataUnique :- grid[j])  )  ))

}

if (stdVar) {
	h = h :* scale_temp
}
return(h)

}
mata mosave lpdensity_bwMSE(), replace
end

********************************************************************************
* imse-dpi bandwidth
********************************************************************************
capture mata mata drop lpdensity_bwIMSE()

mata
real scalar lpdensity_bwIMSE(real colvector data,
							 real colvector grid,
							 real scalar p, real scalar v,
							 string scalar kernel,
							 real colvector Cweights,
							 real colvector Pweights,
							 real scalar massPoints,
							 real scalar stdVar,
							 real scalar regularize,
							 real scalar nLocalMin,
							 real scalar nUniqueMin){

					
ii = order(data, 1)
data = data[ii, 1]
Cweights = Cweights[ii, 1]
Pweights = Pweights[ii, 1]					
n  = length(data)
ng = length(grid)

dataUnique  = lpdensity_unique(data)
freqUnique  = dataUnique[., 2]
indexUnique = dataUnique[., 3]
dataUnique  = dataUnique[., 1]
nUnique     = length(dataUnique)

if (stdVar) {
	center_temp = mean(data)
	scale_temp  = sqrt(variance(data))

	data = (data :- center_temp) :/ scale_temp
	dataUnique = (dataUnique :- center_temp) :/ scale_temp
	grid = (grid :- center_temp) :/ scale_temp
}

if (massPoints) {
	Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique)
} 
else{
	Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights)	
}

weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n
Cweights       = Cweights :/ sum(Cweights) :* n
Pweights       = Pweights :/ sum(Pweights) :* n

weights_normalUnique = runningsum(weights_normal)[indexUnique]
if (nUnique > 1) { 
	weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) 
}

CweightsUnique = runningsum(Cweights)[indexUnique]
if (nUnique > 1) { 
	CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) 
}

PweightsUnique = runningsum(Pweights)[indexUnique]
if (nUnique > 1) { 
	PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) 
}

// obtain preliminary bandwidth for estimating densities
// this is used for constructing preasymptotic matrices
h1  = lpdensity_bwIROT(data, grid, 2,   1,   kernel, Cweights, Pweights, 1, 1, 1, 20+2+1,   20+2+1)

// obtain preliminary bandwidth for estimating F_p+1
// this is used for constructing F_p+1
hp1 = lpdensity_bwIROT(data, grid, p+2, p+1, kernel, Cweights, Pweights, 1, 1, 1, 20+p+2+1, 20+p+2+1)

// obtain preliminary bandwidth for estimating F_p+2
// this is used for constructing F_p+2
hp2 = lpdensity_bwIROT(data, grid, p+3, p+2, kernel, Cweights, Pweights, 1, 1, 1, 20+p+3+1, 20+p+3+1)

dgp_hat = J(ng, 2, .) // Fp+1 and Fp+2 with normalization constants
const_hat = J(ng, 3, .)

for (j=1; j<=ng; j++) {

	// estimate F_p+2
	index_temp = (abs(data :- grid[j]) :<= hp2)
	Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp2 
	Xh_p_temp = J(sum(index_temp), p+4, 1)	
	Xh_p_temp[., 2] = Xh_temp

	for (jj=3; jj<=p+4; jj++) { 
		Xh_p_temp[., jj] = Xh_temp :^ (jj-1)
	}
		
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp)) :/ hp2
	if (kernel == "uniform") 		Kh_temp = 0.5 / hp2
	if (kernel == "epanechnikov") 	Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp2
	
	Kh_temp = select(Pweights, index_temp) :* Kh_temp
	Y_temp = select(Fn, index_temp)
	dgp_hat[j, 2] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+3, 1]  / (hp2^(p+2))

	// estimate F_p+1
    index_temp = (abs(data :- grid[j]) :<= hp1)
	Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp1 
	Xh_p_temp = J(sum(index_temp), p+3, 1)
	Xh_p_temp[., 2] = Xh_temp
	
	for (jj=3; jj<=p+3; jj++) { 
		Xh_p_temp[., jj] = Xh_temp:^(jj-1)
	}
	
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp))  :/ hp1
	if (kernel == "uniform") 		Kh_temp = 0.5 / hp1
	if (kernel == "epanechnikov") 	Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp1
	
	Kh_temp = select(Pweights, index_temp) :* Kh_temp
	Y_temp = select(Fn, index_temp)
	dgp_hat[j, 1] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+2, 1] / (hp1^(p+1))

	// prepare for estimating matrices
    index_temp = (abs(data :- grid[j]) :<= h1)
    Xh_temp = (select(data, index_temp) :- grid[j]) :/ h1 
	
	if (kernel == "triangular")  	Kh_temp = (1 :- abs(Xh_temp)) :/ h1
	if (kernel == "uniform") 		Kh_temp = 0.5 / h1
	if (kernel == "epanechnikov") 	Kh_temp = 0.75 :* (1 :- Xh_temp:^2) :/ h1
	Pweights_temp = select(Pweights, index_temp)

    // estimate Cp matrix
	Xh_p_temp = J(sum(index_temp), (2*p+1) - (p+1) + 1, 1)
	for (jj=1; jj<=(2*p+1)-(p+1)+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(p+jj)
	}
	C_p_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n

	// estimate Cp+1 matrix
	Xh_p_temp = J(sum(index_temp), (2*p+2) - (p+2) + 1, 1)
	for (jj=1; jj<=(2*p+2)-(p+2)+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(p+jj+1)
	}
	C_p1_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n
	
	// estimate S matirx
 	Xh_p_temp = J(sum(index_temp), p - 0 + 1, 1)
	for (jj=1; jj<=p-0+1; jj++) {
		Xh_p_temp[., jj] = Xh_temp:^(jj-1)
	}

	S_hat = cross(Xh_p_temp, Kh_temp :* Pweights_temp, Xh_p_temp) :/ n
	S_hat_inv = invsym(S_hat)
	
	// estimate G matrix
	if (v == 0) {
		G_hat = cross(Xh_p_temp, Kh_temp:^2 :* Pweights_temp, Xh_p_temp) :/ n
	} 
	else {
		if (massPoints) {
			Y_temp = Fn[indexUnique]
			Xh_temp = (dataUnique :- grid[j]) :/ h1 
			Xh_p_temp = J(nUnique, p - 0 + 1, 1)
			Xh_p_Kh_Pweights_temp = Xh_p_temp
		
			if (kernel == "triangular")  	Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp[indexUnique]
			if (kernel == "uniform") 		Kh_temp =  (0.5 / h1) :* index_temp[indexUnique]
			if (kernel == "epanechnikov") 	Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp[indexUnique]
		
			for (jj=1; jj<=p-0+1; jj++) {
				Xh_p_temp[., jj] = Xh_temp:^(jj-1)
				Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique
			}
		
			F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n
		
			G = J(n, p - 0 + 1, .)
		
			for (jj=1; jj<=p-0+1;jj++) {
				G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal
			}

			G_hat = cross(G, G) :/ n
		}
		else {
			Y_temp = Fn
			Xh_temp = (data :- grid[j]) :/ h1 
			Xh_p_temp = J(n, p - 0 + 1, 1)
			Xh_p_Kh_Pweights_temp = Xh_p_temp
		
			if (kernel == "triangular")  	Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp
			if (kernel == "uniform") 		Kh_temp =  (0.5 / h1) :* index_temp
			if (kernel == "epanechnikov") 	Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp
		
			for (jj=1; jj<=p-0+1; jj++) {
				Xh_p_temp[., jj] = Xh_temp:^(jj-1)
				Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights
			}
		
			F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n
		
			G = J(n, p - 0 + 1, .)
		
			for (jj=1; jj<=p-0+1;jj++) {
				G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal
			}

			G_hat = cross(G, G) :/ n
		}	
	}


	// now get all constants
    const_hat[j, 1] = factorial(v) * (S_hat_inv * C_p_hat )[v+1]
    const_hat[j, 2] = factorial(v) * (S_hat_inv * C_p1_hat)[v+1]

    if (v > 0) {
		const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (n * h1))
    } 
	else {
		temp_ii = min((  max((  mean(data :<= grid[j]), 1/n  )), 1 - 1/n ))
		const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (0.5 * n^2) * 
			h1 * temp_ii * (1 - temp_ii))
    }
}

// bandwidth
S = optimize_init()
optimize_init_evaluator(S, &lpdensity_optFunc())
optimize_init_params(S, 1e-8)
optimize_init_argument(S, 1, dgp_hat)
optimize_init_argument(S, 2, const_hat)
optimize_init_argument(S, 3, p)
optimize_init_argument(S, 4, v)
optimize_init_argument(S, 5, "imse-dpi")
optimize_init_which(S, "min")
optimize_init_evaluatortype(S, "gf0")
optimize_init_tracelevel(S, "none")
optimize_init_conv_maxiter(S, 20)
(void) optimize(S)
h = optimize_result_params(S)

if (regularize == 1) {
	for (j=1; j<=ng; j++) {
		if (nLocalMin > 0) {
			h = max((h, sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))]))
		}
		if (nUniqueMin > 0) {
			h = max((h, sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))]))
		}
	}
}

h = min(( h, max((  abs(max(dataUnique)-min(grid)), abs(min(dataUnique)-max(grid))  ))  ))

if (stdVar) {
	h = h :* scale_temp
}
return(h) 

}
mata mosave lpdensity_bwIMSE(), replace
end

********************************************************************************
* optimization criterion function
********************************************************************************
capture mata mata drop lpdensity_optFunc()

mata
void lpdensity_optFunc(todo, h, mat1, mat2, p, v, bwselect, fn, gn, hn) {
	if (bwselect == "imse-dpi") {
		if (v > 0) {
			fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :* mat2[., 1] :+ h :* mat1[., 2] :* mat2[., 2]):^2 + mat2[., 3]:^2 :/ h^(2*v - 1) )
		} else {
			fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :* mat2[., 1] :+ h :* mat1[., 2] :* mat2[., 2]):^2 + mat2[., 3]:^2 :/ h )
		}
	}
	else if (bwselect == "imse-rot") {
		if (v > 0) {
			fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :+ h :* mat1[., 2]):^2 + mat2[., 1]:^2 :/ h^(2*v - 1) )
		} else {
			fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :+ h :* mat1[., 2]):^2 + mat2[., 1]:^2 :/ h )
		}
	} else if (bwselect == "mse-dpi") {
		if (v > 0) {
			fn = h^(2*p+2-2*v) * (mat1[1] * mat2[1] + h * mat1[2] * mat2[2])^2 + mat2[3]^2 / h^(2*v - 1)
		} else {
			fn = h^(2*p+2-2*v) * (mat1[1] * mat2[1] + h * mat1[2] * mat2[2])^2 + mat2[3]^2 / h
		}
	} else {
		if (v > 0) {
			fn = h^(2*p+2-2*v) * (mat1[1] + h * mat1[2])^2 + mat2[1]^2 / h^(2*v - 1)
		} else {
			fn = h^(2*p+2-2*v) * (mat1[1] + h * mat1[2])^2 + mat2[1]^2 / h
		}
	}
	
}
mata mosave lpdensity_optFunc(), replace
end
